/**
 * \copyright
 * Copyright (c) 2012-2017, OpenGeoSys Community (http://www.opengeosys.org)
 *            Distributed under a Modified BSD License.
 *              See accompanying file LICENSE.txt or
 *              http://www.opengeosys.org/project/license
 *
 * \file   WaterDensityIAPWSIF97Region1.h
 *
 * Created on December 8, 2016, 4:19 PM
 */

#pragma once

#include <string>

#include "Material/Fluid/GibbsFreeEnergy/DimensionLessGibbsFreeEnergyRegion1.h"

namespace MaterialLib
{
namespace Fluid
{
/** Water density model
 *  base on the IAPWS Industrial Formulation 1997
 *  <a href="http://www.iapws.org/relguide/IF97-Rev.pdf">IF97-Rev</a>
 */
class WaterDensityIAPWSIF97Region1
{
public:
    WaterDensityIAPWSIF97Region1()
        : _gibbs_free_energy(DimensionLessGibbsFreeEnergyRegion1()),
          _ref_T(1386.0),
          _ref_p(1.653e7),
          _sR(461.526)
    {
    }

    /**  Get density value.
         \param p   Pressure.
         \param T   Temperature.
    */
    double getValue(const double p, const double T)
    {
        const double tau = _ref_T / T;
        const double pi = p / _ref_p;

        return _ref_p / (_sR * T * _gibbs_free_energy.get_dgamma_dpi(tau, pi));
    }

    /**
     *  Get the partial differential of the density with respect to temperature.
     *  \param p   Pressure.
     *  \param T   Temperature.
     */
    double getdValuedT(const double p, const double T)
    {
        const double tau = _ref_T / T;
        const double pi = p / _ref_p;

        const double dgamma_dpi = _gibbs_free_energy.get_dgamma_dpi(tau, pi);
        return -(_ref_p - tau * _ref_p *
                              _gibbs_free_energy.get_dgamma_dtau_dpi(tau, pi) /
                              dgamma_dpi) /
               (_sR * T * T * dgamma_dpi);
    }

    /**
     *  Get the partial differential of the density with respect to pressure.
     *  \param p   Pressure.
     *  \param T   Temperature.
     */
    double getdValuedp(const double p, const double T)
    {
        const double tau = _ref_T / T;
        const double pi = p / _ref_p;

        const double dgamma_dpi = _gibbs_free_energy.get_dgamma_dpi(tau, pi);
        return -_gibbs_free_energy.get_dgamma_dpi_dpi(tau, pi) /
               (_sR * T * dgamma_dpi * dgamma_dpi);
    }

private:
    const DimensionLessGibbsFreeEnergyRegion1 _gibbs_free_energy;

    const double _ref_T;  ///< reference temperature in K.
    const double _ref_p;  ///< reference pressure in Pa.
    /// Specific water vapour gas constant in J/(kgK).
    const double _sR;
};

}  // namespace Fluid
}  // namespace MaterialLib
